1 Showcasing RRPlots

1.0.1 Libraries

library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
layout(matrix(1:1, nrow=1))

1.0.2 Wisconsin Data Set

dataBreast <- read.csv("~/GitHub/RISKPLOTS/DATA/wpbc.data", header=FALSE)
table(dataBreast$V2)
## 
##   N   R 
## 151  47
rownames(dataBreast) <- dataBreast$V1
dataBreast$V1 <- NULL
dataBreast$status <- 1*(dataBreast$V2=="R")
dataBreast$V2 <- NULL
dataBreast$time <- dataBreast$V3
dataBreast$V3 <- NULL
dataBreast <- sapply(dataBreast,as.numeric)
## Warning in lapply(X = X, FUN = FUN, ...): NAs introduced by coercion
dataBreast <- as.data.frame(dataBreast[complete.cases(dataBreast),])
table(dataBreast$status)
## 
##   0   1 
## 148  46

1.1 Exploring Raw Features with RRPlot

convar <- colnames(dataBreast)[lapply(apply(dataBreast,2,unique),length) > 10]
convar <- convar[convar != "time"]
topvar <- univariate_BinEnsemble(dataBreast[,c("status",convar)],"status")
pander::pander(topvar)
V35 V24 V34 V7 V16 V14 V17
0.0261 0.0261 0.0261 0.0623 0.126 0.126 0.126 
topv <- min(5,length(topvar))
topFive <- names(topvar)[1:topv]
RRanalysis <- list();
idx <- 1
for (topf in topFive)
{
  RRanalysis[[idx]] <- RRPlot(cbind(dataBreast$status,dataBreast[,topf]),
                              atProb=c(0.90),
                  timetoEvent=dataBreast$time,
                  title=topf,
#                  plotRR=FALSE
                  )
  idx <- idx + 1
}

names(RRanalysis) <- topFive

1.2 Reporting the Metrics

ROCAUC <- NULL
CstatCI <- NULL
RRatios <- NULL
LogRangp <- NULL
Sensitivity <- NULL
Specificity <- NULL

for (topf in topFive)
{
  CstatCI <- rbind(CstatCI,RRanalysis[[topf]]$c.index$cstatCI)
  RRatios <- rbind(RRatios,RRanalysis[[topf]]$RR_atP)
  LogRangp <- rbind(LogRangp,RRanalysis[[topf]]$surdif$pvalue)
  Sensitivity <- rbind(Sensitivity,RRanalysis[[topf]]$ROCAnalysis$sensitivity)
  Specificity <- rbind(Specificity,RRanalysis[[topf]]$ROCAnalysis$specificity)
  ROCAUC <- rbind(ROCAUC,RRanalysis[[topf]]$ROCAnalysis$aucs)
}
rownames(CstatCI) <- topFive
rownames(RRatios) <- topFive
rownames(LogRangp) <- topFive
rownames(Sensitivity) <- topFive
rownames(Specificity) <- topFive
rownames(ROCAUC) <- topFive

pander::pander(ROCAUC)
  est lower upper
V35 0.650 0.555 0.744
V24 0.634 0.542 0.725
V34 0.661 0.577 0.746
V7 0.610 0.515 0.705
V16 0.598 0.504 0.692 
pander::pander(CstatCI)
  mean.C Index median lower upper
V35 0.622 0.622 0.518 0.712
V24 0.677 0.676 0.595 0.751
V34 0.659 0.659 0.591 0.723
V7 0.666 0.663 0.583 0.742
V16 0.614 0.616 0.525 0.698 
pander::pander(RRatios)
  est lower upper
V35 1.37 0.698 2.69
V24 1.93 1.122 3.31
V34 1.44 0.741 2.82
V7 1.37 0.700 2.69
V16 1.04 0.462 2.32 
pander::pander(LogRangp)
V35 0.16213
V24 0.00331
V34 0.37299
V7 0.03818
V16 0.50639 
pander::pander(Sensitivity)
  est lower upper
V35 0.152 0.0634 0.289
V24 0.239 0.1259 0.388
V34 0.152 0.0634 0.289
V7 0.152 0.0634 0.289
V16 0.109 0.0362 0.236 
pander::pander(Specificity)
  est lower upper
V35 0.899 0.838 0.942
V24 0.899 0.838 0.942
V34 0.899 0.838 0.942
V7 0.899 0.838 0.942
V16 0.899 0.838 0.942 
meanMatrix <- cbind(ROCAUC[,1],CstatCI[,1],Sensitivity[,1],Specificity[,1],RRatios[,1])
colnames(meanMatrix) <- c("ROCAUC","C-Stat","Sen","Spe","RR")
pander::pander(meanMatrix)
  ROCAUC C-Stat Sen Spe RR
V35 0.650 0.622 0.152 0.899 1.37
V24 0.634 0.677 0.239 0.899 1.93
V34 0.661 0.659 0.152 0.899 1.44
V7 0.610 0.666 0.152 0.899 1.37
V16 0.598 0.614 0.109 0.899 1.04

1.3 Modeling

ml <- BSWiMS.model(Surv(time,status)~1,data=dataBreast,NumberofRepeats = 10)

[++++++++++++++++++++++++++++++++++++++++++++++++++++++]…..

sm <- summary(ml)
pander::pander(sm$coefficients)
Table continues below
  Estimate lower HR upper u.Accuracy r.Accuracy
V24 4.49e-02 1.01 1.05 1.08 0.598 0.263
V26 5.33e-03 1.00 1.01 1.01 0.593 0.237
V27 2.36e-04 1.00 1.00 1.00 0.608 0.273
V34 1.14e-02 1.00 1.01 1.02 0.634 0.259
V7 5.75e-08 1.00 1.00 1.00 0.588 0.313
V35 2.13e-03 1.00 1.00 1.00 0.727 0.603
V6 9.72e-08 1.00 1.00 1.00 0.577 0.237
V4 1.88e-06 1.00 1.00 1.00 0.582 0.237
Table continues below
  full.Accuracy u.AUC r.AUC full.AUC IDI NRI z.IDI
V24 0.598 0.609 0.508 0.609 0.0618 0.435 2.86
V26 0.593 0.598 0.500 0.598 0.0626 0.393 2.77
V27 0.609 0.608 0.510 0.608 0.0563 0.435 2.76
V34 0.632 0.618 0.507 0.617 0.0315 0.468 2.40
V7 0.591 0.595 0.522 0.599 0.0475 0.378 2.26
V35 0.607 0.641 0.609 0.602 0.0284 0.551 2.26
V6 0.577 0.588 0.500 0.588 0.0459 0.353 2.19
V4 0.582 0.606 0.500 0.606 0.0449 0.426 2.19
  z.NRI Delta.AUC Frequency
V24 2.65 0.10154 1.0
V26 2.38 0.09827 1.0
V27 2.64 0.09730 1.0
V34 2.83 0.11073 1.0
V7 2.29 0.07623 0.9
V35 3.41 -0.00641 0.8
V6 2.13 0.08813 0.1
V4 2.62 0.10649 0.2

1.4 Cox Model Performance

Here we evaluate the model using the RRPlot() function.

1.4.1 The evaluation of the raw Cox model with RRPlot()

Here we will use the predicted event probability assuming a baseline hazard for events withing 5 years

index <- predict(ml,dataBreast)
timeinterval <- 2*mean(subset(dataBreast,status==1)$time)

h0 <- sum(dataBreast$status & dataBreast$time <= timeinterval)
h0 <- h0/sum((dataBreast$time > timeinterval) | (dataBreast$status==1))
pander::pander(t(c(h0=h0,timeinterval=timeinterval)),caption="Initial Parameters")
Initial Parameters
h0 timeinterval
0.323 51.1 
rdata <- cbind(dataBreast$status,ppoisGzero(index,h0))
rownames(rdata) <- rownames(dataBreast)

rrAnalysisTrain <- RRPlot(rdata,atProb=c(0.90,0.80),
                     timetoEvent=dataBreast$time,
                     title="Raw Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

1.4.2 Uncalibrated Performance Report

pander::pander(t(rrAnalysisTrain$OERatio),caption="O/E Ratio")
O/E Ratio
est lower upper
0.844 0.618 1.13 
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
1.02 1.02 0.964 1.07 
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
0.797 0.797 0.789 0.806 
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
mean.C Index median lower upper
0.682 0.681 0.597 0.754 
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.64 0.548 0.731 
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.261 0.143 0.411 
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.838 0.942 
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.418 0.361 
pander::pander(t(rrAnalysisTrain$RR_atP),caption="Risk Ratio")
Risk Ratio
est lower upper
1.84 1.08 3.14 
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 12.346960 on 2 degrees of freedom, p = 0.002084
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 147 29 36.99 1.725 8.974
class=1 20 5 4.11 0.193 0.216
class=2 27 12 4.90 10.269 11.609